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ABSTRACT 

This article synthesizes observational data from an extensive program aimed toward a comprehensive 
understanding of star formation in a low-mass star-forming molecular cloud. New observations and 
published data spanning from the centimeter wave band to the near infrared reveal the high and low 
density molecular gas, dust, and pre-main sequence stars in L1551. 

The total cloud mass of ~ 160 contained within a 0.9 pc has a dynamical timescale, idyn = 
1.1 Myr. Thirty-five pre-main sequence stars with masses from ~ 0.1 to 1.5 Mq are selected to be 
members of the L1551 association constituting a total of 22 ± 5 Mq of stellar mass. The observed star 
formation efficiency, SFE = 12%, while the total efficiency, SFEtot, is estimated to fall between 9 and 
15%. 

L1551 appears to have been forming stars for several idyn with the rate of star formation increas- 
ing with time. Star formation has likely progressed from east to west, and there is clear evidence 
that another star or stellar system will form in the high column density region to the northwest of 
L1551 IRS5. 

High-resolution, wide-field maps of L1551 in CO isotopologue emission display the structure of the 
molecular cloud at 1600 AU physical resolution. The ^'^CO emission clearly reveals the disruption of 
the ambient cloud by outfiows in the line core and traces the interface between regions of outflow and 
quiescent gas in the line wings. Kinetic energy from outflows is being deposited back into the cloud 
on a physical scale Apeak ~ 0.05 pc at a rate, i?input ~ 0.05^0. The remaining energy afforded by the 
full mechanical luminosity of outflow in L1551 destroys the cloud or is otherwise lost to the greater 
interstellar medium. 

The C^^O emission is optically thin and traces well the turbulent velocity structure of the cloud. The 
total turbulent energy is close to what is expected from virial equilibrium. The turbulent velocities 
exist primarily on small scales in the cloud and the energy spectrum of turbulent fluctuations, E(k) oc 
fc"'^, is derived by various methods to have (3 ~ 1-2. The turbulent dissipation rate estimated using 
the results of current numerical simulations is -Ediss ~ ^-input- 

This study reveals that stellar feedback is a significant factor in the evolution of the L1551 cloud. 
Subject headings: stars: formation — ISM: clouds — stars: pre-main sequence — ISM: individ- 
ual(L1551) — ISM: structure — radio lines: ISM — techniques: interferometric 



1. INTRODUCTION 

Low-mass stars form from dense, quiesce nt cores em- 
bedded within turbulent molecular clouds ( Myers et al. 
[1983, Beichman et al.||1986| lLarson|p8l] ). Followiiigi :' 
rapid inside-out gravitational coliapse phase ( Shu||1977 1 , 
embedded sources accrete material from their surround- 
ings through a di sk while bipolar flows interact with the 
progenitor cloud (Harti gan et al.|1995 Richer et al. 2000 



4:Mq yr ^ (Prantzos & Aubert 19951, indicates that 



[Reipu rth & Bally 2001 Young stars become optically 
visible and accretion diminishes over time, as does the 



mass of the circumstellar disk (Hartmann et al. 1998 



Lada 1991 1. The luminosities of the pre-main sequence 



stars evolve as they contract and e ventually settle onto 
the main sequence (Hayashi 19611. This succession in 



the formation of low-mass stars is fairly certain, but the 
theoretical context in which these stages fit is still de- 
bated. 

As first pointed out by Zuckerman fc Evans ( 1974[), 
the m olecula r m ass in the Galaxy, ~ 2 x W^Mq ( Bloe- 
men et aL]|l986[ ), compared to the star formation rate. 



star formation is a glob ally inefficient process (see also 
Krumholz fc: Tan 2007[ ). Another universal property of 
star formation is the distribution of stellar masses pro- 
duced by a star-forming cloud, or the initial mass func- 
tion (IMF; [Meyer et "aLl|2000i. Any successful theory 
of star formation will explain these general properties of 
star formation in a way consistent with the above sce- 
nario. 

In th e standard paradigm for low-mass star formation 
(Shu et al. 1987j), magnetically sub-critical cores evolve 
into pr otostars ov er timescales of ~ 10 Myr (Ciolek & 



MouschoviaS|| 1994[ ) . Upon the formation of a protostar 
outflows limit the amount of material available for accre- 
tion and disperse the parent cloud. Thus feedback from 
young sources play a critical role in determining the flnal 
stellar masses and overall star formation efficiency. The 



copious energies of molecular outflows (.Bachiller 1996 1, 
the large sp atial extent of Herbig-Haro flows ( [Reipurth &z 
Bally 20011) , and the existence of remnant cavities m star- 
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formi ng clouds (Quillen et al. 2005 Cunningham et al. 
20061 all lend credence to this picture. However, mag- 



netic field observations do not show strong evi dence for 
magnetically sub-critical cores ( |Crutcher|[l999l ). 

In another description of low-mass star formation, or- 
dered magnetic fields and stellar feedback are not critical 
elements of the process. Rather star formation occurs on 
the order of a dynamical time scale in regions of converg 



ing turbulent velocity fields ( Ballesteros-Paredes et al. 
2007 lElmegreen 20001. The inetticiency oi star torma- 



tion and the IMF arise solely from the properties of inter 



stellar turbulence in this framework ( Vazquez-Semadeni 
eTaljpOOS] IPadoan et allpOOT] [Padoan fc JNordlund 



2002[ |. Turbulence in star forming regions is difiicult to 

charac terize from observational data alone fElmegreen' 
& Scalo 2004) and support for this theory relics heavily 
upon computer simulation. 

The aim of this work is to gain insight into the rela- 
tive importance of different mechanisms involved in the 
star formation process by doing a detailed case study 
of a star-form ing region. The Lynds dark cloud L1551 
(Lynds 19621 is an ideal subject for this kind of study 
since it IS a relatively isolated, nearby 160 pc), star- 
forming cloud with a wealth of activity that encom- 
passes all phenomena known to be associated with low- 
mass star forma tion includin g; a population of pre-main 
sequence stars (iKenyon fc^ Hartmann 1995), multiple 
outflow s ('Snelll'ra ^ |Moriarty-Schieven fc S ncU 1988; 
Moriarty-Schicvcn et al.||1995| [Pound fc Bally 1991) 



jets ([M undt et al. 1990*), wmasf Welc h et aTpD OO Giova 

iogas 



nardi et al. 2000), an abundance of shocked gas (Dcvine 
et al. 1999]), and reflection nebulosity (Draper ct al. 
T^85) . This article is based on the work by ,Swift, ( ,2006 ) , 
and summarizes some of the conclusions therein while 
elaborating on some of the themes via new data and anal- 
yses. 

Our new observations are presented in §[2] and con- 
tribute significantly to the numerous observations in this 
region. Data pertaining to the young stars around L1551 
give insights into the star formation history and the cloud 
lifetime, and Spitzer IRAC observations of the dense core 
L1551-MC add to the understanding of the highest col- 
umn regions of L1551 in §l3] The general properties of 
the LI 551 cloud are then deduced in §|4] using dust ex- 
tinction and ^■^CO emission. 

The high-resolution CO isotopologue maps are shown 
in §|5]whcre the ^'^CO emission proves to be a good tracer 
of stellar feedback in LI 551. The nature of this feedback 
is explored in §[6] where it is found that outfiow in L1551 
is both stirring the ambient gas as well as excavating 
mass from the cloud. The C^^O emission traces the mass 
distribution in L1551 better than the ^'^CO, and these 
data are used to probe the turbulent velocity field in 
the cloud in §J7] The star formation history in L1551 is 
discussed in §[8[ and the main conclusions of the article 
are summarized in §[9j 

2. OBSERVATIONS AND REDUCTIONS 

2.1. Wide-Field BIMA Mosaic 

A 250 pointing intcrfcrometric mosaic conducted with 
the BIMA^ interferometer covers ~ 140 square arcmin- 

^ The Berkeley-IUinois-Maryland Association array was, at the 
time of these observations, operated by the University of Cahfornia, 



6 dWelch et 
windows 



utes of the L1551 cloud. The mosaic is designed to image 
the regions of highest ^"^CO emis sion as seen in the single- 
dish maps of Swift et al. ( 2005 1 including the region of 
L1551-MC. Tne full mosaic is broken into 13 sub-mosaics, 
labeled A-M, each consisting of roughly 20 pointings. 
Each sub- mosaic track was run for a full night in both the 
C and D array configurations while follow-up scheduling 
was directed to create nearly constant sensitivity across 
the map. Table 5.1 in Swift (20061 summarizes the ob- 
servations which span 7 years with the vast majority of 
tracks run between 2001 and 2004. Figure [T| a) displays 
the pointing centers for the s ub-mosaics superp osed on 
velocity integrated "CO from [Swift et al.| ( |2005| . 

The pointing coordinates in each sub-mosaic are de- 
fined in terms of right ascension and declination offsets 
from the most central pointing in the sub-mosaic. The 
offsets are chosen such that the pointing centers make a 
grid of equilateral triangles with each side equal to the 
Nyquist sampling angle, Gn = 46". The exception is 
Mosaic M where both the right ascension and declina- 
tion offsets are integer multiples of On. 
The mosaic data w ere taken in correlator mode 
al. 1996 1 with the two high-resolution 
containing 256 channels set to 12.5 MHz 
widths centered on the transition frequencies of C^^O 
(109.78216 GHz) and ^^CO (110.20135 GHz) in the up- 
per side band. This translates to a velocity resolution of 
0.13 km s^^ in these spectral windows. 

Each sub-mosaic track cycled through all pointings 
with 23 s integrations performed in groups of 1, 2 or 3 
at each pointing center depending on mosaic size in con- 
sideration of proper phase calibration. Integrations at 
all pointings were performed between phase calibration 
for tracks with 20 pointings or less. For the larger mo- 
saics (A, B, and L), two calibration integrations were 
embedded in the mosaic to ensure phase tracking. In 
our final data, the u—v sampling is complete over the en- 
tire map from 2 < + ^ 25 kA. The observations 
were phase calibrated using nearby quasars 0530-1-135 
and 0429+113, both > IJy sources at these observing 
frequencies throughout the program. 

The flux scale was calibrated using W3(0H) as the 
primary calibrator. Saturn or Uranus was used as a sec- 
ondary calibrator, and Jupiter was used occasionally for 
D configuration observations. The fiux scale in the final 
maps is estimated to be accurate at the ~ 15% level. 

Anomalous amplitudes and phases were fiagged from 
the data, and no shadowed data were used. Line- 
length and gain calibrations were applied, and then 
all data were Fourier trans formed simultaneously using 
MIRIAD ( [Sault et al.|1995[ ). Natural weighting was used 
to maximize signal to noise. Figure ITT &) shows the pos- 
itive values of the velocity integratea emission from the 
Fourier transform — the "dirty" map. 

The maps were then deconvol ved with a Steer- 
Dewdney-Ito CLEAN algorithm (St eer et al.|1984 
a low loop gain of 0.05 and a maximal number o: 
tions so that all fiux above the 1 a level was transferred to 
CLEAN components. The clean image was then restored 



using 
itera- 



Berkeley, the University of Ilhnois, and the University of Maryland 
with support from the National Science Foundation. The BIMA 
array has now been combined with the Owens Valley Radio Obser- 
vatory's millimeter interferometer to create the Combined Array 
for Research in Millimeter- Wave Astronomy, or CARMA. 
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with a bivariate Gaussian clean beam having FWHM di- 
mensions of 10."3 X 8."7 and position angle of 29°. 

2.2. Combining Single-Dish and Interferometer Data 
A single pointing of an interferometer cannot measure 
visibilities where \/u'^ + v'^ < D. Observations using 
multiple pointings allow the u-v plane to be sampled 
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Fig. 1.— (a) Velocity integrated l^CO(l - 0) emission from 
the Kitt Peak 12 m showing the pointings from each sub-mosaic 
A-M. (b) The integrated emission from th e int erferometer data, 
and (c) the combined map as described in § |2.2| The gray scale is 
transferred linearly between and 15 K km s^^ for all images. 



slightly inside this limit due to the Ekers & Rots ( 1979 1 
effect, but intensity variations on scales larger than the 
primary beam cannot be recovered with interferometry 
alone [compare Figure [ij a) and (&)]. Therefore, when 
mapping extended objects with interferometry the addi- 
tion of single dish data is essential in reconstructing ac- 
curate depiction s of the sky brightness distribution (see 
Holdaway||1999| . 

Approximately 6% of the total ^"^CO flux is recovered 
with our cleaned BIMA mosaic. The flux recovery is 
greatest in the line wings, reaching as high as ~ 35%, 
and steadily decreasing to effectively zero in the line core. 
The flux recovery in the C^^O map is on the level of a 
few percent. 

The flux density scales of the single-dish and interfero- 
metric datasets were compare d in the region of u-v over - 
lap in the manner described by Stanimirovic et al. ( 1999 1. 
The flux density scale for the Kitt Peak 12 m was con- 
sisten t with the nominal scale of 33 Jy (Heifer et al. 
20031, and a flux scaling factor of unity is assumed tor 
all data in the combination. 

The data are combined using a linear combination of 
the two d atasets followed by a maximum entro py decon- 
volution ( Stanimirovic et ar]|1999 Holdaway|1 999 ) . The 
regridded single-dish map, seen in velocity integrated 
^■^CO intensity in Figure |l|a), is flrst combined with 
the dirty BI MA mosaic using Equation 4 in |Stanimirovic| 
et al. (19991. The ratio of beam areas, a = 0.029. 

TEe maximum entropy deconvolution was performed 
using MIRIAD. An rms factor of 1.2 was used to ac- 
count for the discrepancy between the theoretical and 
the real rms in the images. The regridded single-dish 
image served as the default image, and the total flux of 
the output was forced to agree with the total flux in the 
default image. Without using the single-dish image as 
the default the total flux in the final image was overesti- 
mated by nearly a factor of 2. 

Convergence was achieved in all 64 image planes in < 
15 iterations. There are no significant artifacts from the 
deconvolution and the overall image quality is very good. 
The deconvolution residuals contain minimal structure 
with a consistent rms from channel to channel, evidence 
of a successful deconvolution. The regridded single-dish 
map, the interferometer map, and the deconvolved com- 
bined map integrated over the ^"^CO line profile can be 
seen in Figure [T] 

2.3. Deep Near Infrared Imaging and A Large Scale 
Extinction Map 

Deep near infrared (NIR) imaging was obtained on 
the night of 2004 November 11 to s upplement data 
from the 2MASS point source catalog ( Skrutskie et al. 
20061 in the L1551 region. The images were obtained 
with FLAMINGOS on the Kitt Peak 2.1m telescope. 
FLAMINGOS is a combination wide-field, near infrared 
imager and multi- object spectromete r built by the Uni- 
versity of Florida (Elston et al. 2003). In imaging mode, 
the 2048 x 2048 science grade Hawaii-II HgCdTe array 
has a plate scale of 0."606 pixel" ^ and a 20.'7 field of view. 

Three sets of exposures were taken through the J, H, 
and X-band filters with total exposure times of 720, 720, 
and 875 s consisting of 15 45 s exposures in J and H, and 
25 35 s exposures in K. These data were red uced using 
the LONGLEGS (|Roman-Zuniga et al.|r2006[) reduction 
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procedure. 

J, H , and band photometric data from the 2MASS 
point source catalog were selected within a 4 ° radius of 
LI 551 using VizieR (Ochsenbein et al. 20001. Aperture 



photometry of 171 FLAMINGOS point sources matched 
with the 2MASS catalog produced zero points for the 
FLAMINGOS data to an accuracy of - 0^25. The hm- 
iting magnitudes of our FLAMINGOS data are 22.3 in 

J, 21.9 in H, and 21.4 in K^. 

Alves '2001|) is an 



The NICER method (Lombardi 



optimized technique for converting NIR color excess of 
background stars into an Ay map. Using conversion 



factors from the literature (Harvey et al. 2003 Rieke & 



Lebofsky 1985), extinction values were assigned to grid- 
ded 1.'2 pixels based on a weighted mean of measured ex- 
tinctions from stars falling within a 2. '9 FWHM Gaussian 
beam centered on the pixel coordinates. The weighting 
incorporates errors in photometry, errors due to intrinsic 
variations of stellar colors determined from a control field 
located at 4^20^00% 17°00'00" (J2000), and the stellar 
position relative to the pixel center. 

Figure |2] shows the extinction map of L1551 with a 
dynamic range of 0"!5 < Ay < 18™ and a mean (Ay) « 6 
taken over the central 20' x 20' of the cloud. The rms 
in the map varies from 0™6 in the part of the map with 
only 2MASS data, to 0™3 in the outskirts of the Kitt 
Peak 2.1m data, to more than 1™ in the inner regions 
where there are a paucity of background stars. The errors 
are skewed toward higher extinction values in regions of 
high obscuration. 



2.4. Spitzer IRAC Imaging 
A single I nfrared Array Camera (IRAC; 



Fazio 



et al. 20041 footprint centered on the position 
of peak NH3 emiss ion from L1551 -MC, 4'^31™09".9, 



+18°12'41" (J2000) ( |Swift et al.|2005| , was imaged with 
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Fig. 2. — Map of the measured dust extinction from tlie L1551 
cloud in units of Ay ■ Known pre-main sequence stars are overlaid 
for reference, and the 2'.9 Gaussian smoothing kernel is shown at 
the bottom left. 



the Spitzer Space Telescope on 23 March 2006. In 8min 
of total on source integration time, a 3 cr point source 
sensitivity of roughly 2, 4, 10, and 10/iJy was achieved 
in the 3.6, 4.5, 5.8 and 8.0 /zm bands, respectively. The 
post-BCD^ data were inspected and found to be free 
from artifacts in the first three bands. Small, but no- 
ticeable (~ 1 cr peak) decrements in the floor value of 
several north-south columns produce the "jail bar effect" 
in our 8 /im image. Aperture photometry provided point 
source flux measurements (DNs^^) using an aperture ra- 
dius of 4 pixels, or 4'f8. Zeropoint values of 280.9, 179.9, 
115.0, and 64.13 Jy in the 4 respective IRAC bands con- 
vert these fluxes into magnitudes. 

3. THE STELLAR COMPONENT 

Figure [3] shows the velocity integrated CO (1 — 0) emis- 
sion f rom the southern e nd of the Taurus molecular com- 
plex (Dame et al. 20011. To the north, the majority of 
young star s in T^irus are being born along high density 



filaments (Onishi et al. 1998| 



Hartmann 2002}. L1551 
17! 



hes further out of the Galactic plane at I « 179°, and 
appears as a n isolated molecular clump at a distance of 



160 ± 20 pc (Bertout et al. 1999 Snell 



19811 



There is an association of pre-mam sequence (PMS) 
stars centered near the peak CO column density of 
L1551. Table [1] displays 70 PMS stars within a 4° ra- 
dius of L1551. There have been numerous sear c hes for 
PMS stars in this regio n fHerbig fc Belli |1995| [Feigel- 

Briceiio e^al.||199^ 



son et al. 19871 



Gomez et al. 1992 

2002| ILuhniiaiif |2000p and it is expected that all ' PMS 
stars within the central square degree have been identi- 
fied. Many searches for companion stars have also been 
conducted in this region sensitive to separations between 
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Fig. 3. — Southern part of the Taurus molecular complex in 
CO(l — 0) emission shown in gray scale (0-25K kms~^) and con- 
tours (6, 12 and 18Kkms~^) jPame et al.|20o"l] | w ith th e positi ons 
of pre -main sequence stars {star symbols) overlaid | |Pall a & Stali ler| 
|2002| also see Tablem. L1551 is seen as a relatively isolated molec- 
ular clump coincident with an association of young stars. The posi- 
tions of well-known star-forming clouds TMC 1, TMC 2, and L1536 
as well as T Tau are labeled. 



see [http : //ssc . spitzer . caltech.edu/postbcd 
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1 and 2000 AU (Kohler & Leinert 


1998 


Simon et al.|1995 


Ghez et al. 


|1993| ILeinert et al.| 




IReipurth & Zin-| 


necker|1993 


1, and a large fraction of companion stars are 



3.1. Selection 

Stars with a high probability of having formed from the 
LI 551 cloud are chosen from the list in Table [T| accord- 
ing to their proper motions and radial velocities where 
possible, as well as the projected spatial distribution. 

The 34 proper motion entries listed in Table [T] have 
typic al measurement errors of 0.5 arcseconds per cen- 
tury ( Ducourant et al.] |2005) translating to a '--^ 5kms~^ 
error on the two-dimensional space motions. The 
mean proper motion for this sample is (Aa, AS) = 
(0.71±0.1, — 1.5±0.1) arcseconds per century. There are 
three outliers in the distribution of proper motions that 
lie beyond 15kms-i of the mean— RX J0433.7-I-1823, 
RXJ0430.8-I-2113, and RX J0444.4-H1952— and are re- 
jected from the sample. 

The distribution of 21 radial velocities have a mean 
of 5kms^^ (LSR) and a dispersion consistent with the 
quoted measurement errors for this sample, 2kms~^ 
plerbig & Bell 1995). The only clear outlier in this dis- 
tribution is DQ Tau, and it is rejected based on this cri- 
terion. 

Figure [4] shows the projected radial distribution of 
all PMS stars from Table [T] The stellar density drops 
smoothly to zero and then persists at a roughly constant, 
low level beyond 2°. Based on this distribution, all stars 
beyond a 1?6 radius from LI 551 are rejected from our 
sample. 

Of the spatially selected PMS population, 5 have un- 
known spectral types and 1 has inadequate photometry. 
These stars are flagged from the s ample. Main- sequence 
colors (Kenyon & Hartmann 1995 1 are assumed for weak- 
line T - l'auri stars (wT'T' Ss) a,nd the NIR color-color locus 
from (IMeyer et al.|1997 ) is assumed for classical T-Tauri 



stars (cTT'Ss). Seven stars in the sample have NIR colors 
that are suspect or inconsistent with these assumptions 
and they are marked with a "C" in Table [T] Six of these 
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Fig. 4. — Histogram of the number of pre-main sequence (PMS) 
stars within 0?25 bins of projected radius from L1551. There is 
an absence of PMS stars between 1?6-2?1, and this gap is used 
to separate the L1551 candidates {cross-hatched area) from PMS 
stars in the field that are less likely to have formed in L1551. 
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Fig. 5. — Hertzsprung-Russell diagram for the 30 stars that 
passed our selection criteria of § |3.1| (embedded sources are ex- 
cluded). The assu med error in lummosity is shown at the bot- 
tom left (see § |3.2| l. Pre-main sequence evolutionary tracks for 
mass es between U.1-1.2M(7) and ages between 1-20 Myr are over- 
laid i jBaraffe et a l. 1998, 2002). The two stars to the upper left, 
TAP 51 and HBO 407, lie near the main sequence {dashed line) and 
are likely interlopers. 



seven have adequate supplementary data from the litera- 
ture that are used in further analyses, while the remain- 
ing source, RX J0437.4-f 1851 A, is flagged. The number 
of stars used in the following analyses totals 37 (39 if 
the binary companions to L1551 IRS5 and L1551 NE 
are counted), 7 (9) of which are embedded. There are 
29 stellar systems and 10 companions in the sample for 
a companion star fraction of 34%. This is lower than, 
but comparable to the measured companion star fraction 



in lar ge samples of young stars {e.g., Kohler &: Leinert 
T998|. 



3.2. The Hertzsprung-Russell Diagram 

The spectral types of the PMS stars are converted to 
effective temperatures using Kenyon & Hartmann' (19951 
Table A5) fo r spectral types earlier than MO and Luhman] 



et al. (2003 Table 8) for later spectral types. Ett'ective" 
temperatures for fractional spectral types are computed 
by interpolation. 

The extinction toward each source is derived by tracing 
the reddening line back to either the cTTS locus or the 
main sequence, otherwise literature values are used. The 
NIR fluxes were then correcte d according to the mod- 
els of Rieke & Lebofsky (19851. Bolometric luminosities 
are computed by applying a bolometric correction to the 
extinction corrected J band magnitude, Jc. Ke nyon "fc] 
Hartmann|(1995, Table A5) is used for spectral types ear- 
ier than MO, Bessell. (1991) Table 2) for spectral types 



MO through Isil and [Rad et al.| ( |2001[ Table 2) for spec- 
tral types M8 through M9.5. 

Figure [5] displays the selected sample of PMS 
stars, excluding embedded sources, on the theoretical 
Hertzsprung-Russel l diagram. A subset of the PMS evo- 



lutionary tracks by Baraffe et al. (1998 2002) are over 



laid. The errors o n the data points in this plot are based 



on arguments by Hartmann (20011. Our assumed error 
of S log L = 0.16 is somewhat conservative since this sug- 
gested value was derived assuming all stars to have unre- 
solved companions while it is expected that most, if not 
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all, companion stars in our final sample are resolved. No 
error is assumed on the effective temperatures since, in 
conjunction with the bolometric correction, these errors 
will have a small effect on the derived ages. 

3.2.1. Derived Stellar Properties 

Table [2] displays the ages and masses of the L15 51 pop- 
ulatio n derived from the evolutionary model s of IBaraffe 



et al. 1(1998 2002) along with the type, class (11 
Adams et arp9 87), and the derived Ay ana 



I Lada|1987 



Jc values. 

i'he errors in the determined age are calculated by re- 
deriving these quantities for logL ± i51ogL. The errors 
in mass are estimated to be < 25% for a half spectral type 
uncertainty. Two stars fall very near the zero-age main 
sequence — HBC 407, and TAP 51 — and are heretofore re- 
jected based on the likelihood that they are interloping 
main-sequence stars. 

Figure [6] shows the cumulative age distribution of pre- 
main sequence stars in L1551. Each bin value represents 
the number of stars associated with L1551 that have ages 
greater than or equal to the age value of the bin, and 
the errors are determined using a Monte-Carlo technique. 
There are a total of 35 stars in this sample, 28 have the- 
oretical ages from Figure [5] and the 7 embedded sources 
assumed to have ages < IMyr fall within the first bin. 
The total stellar mass of this sample is 22 ± 5 AIq . 

The majority of the selected PMS stars have ages less 
than 4Myr. However, both age distributions in Fig- 
ure |6] derived from independent PMS evolutionary mod- 
els show that ^ 20% of the stars in this region formed 
more than 6Myr ago. While the derived age spread in 
this region is larger than the age spread typ ically stated 
for the Taurus complex as a whole (e.g ., Ballesteros- 
Paredes et al.||1999 Hartmann et al. | 2001|, tnese results 



are consistent with previous studies ( [Gomez et al.||1992 
[Palla fc Stahler.,2002^ . 

3.3. The Spatial Distribution of the PMS Population 

The radial distribution of PMS stars in L1551 is con- 
centrated near the center of the cloud, and decreases to 
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Fig. 6. — Cumulative histogram of pre-main sequence ages for the 
L1551 association derived from the Hertzsprung-Russell diagram. 
Each 2 Myr bin represents the number of stars with ages greater 
than or equal to the bin age. The solid line represents the ages 
obtained with models of Baraffe ct al. ( 1998 , 2002 1 and the dotted 
line represents the ages obtained with the , Palla fc Stahler| | )1999| ) 
models. 



zero at a radius of -Rstais = 1°6, or 4.5 pc in projected 
distance. The dispersion timescale for the stellar associ- 
ation is estimated 

,^ -^stars ^-|^^ 



disp 



'3D 



where Ugj-, is the three-dimensional velocity dispersion of 
the stars, and the factor of tt/2 accounts for a random 
inclination angle of stellar motion with respect to the 
line of sight. The dispersions of the proper motion and 
radial velocity distributions of the L1551 association are 
approximately equal to the stated m easurement errors 
( [Ducourant et al.||2005[ [Herbig fc Bell,, 1995 j . Therefore 
an upper limit on the one-dimensional velocity dispersion 
for the stellar association of 2kms~^ results from the 
dispersion of radial velocities. A lower limit of 0.3 km 
is derived from the velocity dispersion of C^^O emission 
m the cloud (§[7|. 
For a one-dimensional velocity dispersion of Ikms^^ 



disp 



4 Myr. The observed spatial 
distribution of stars is thus consistent with dispersion 
by random motions over several million years. It is also 
possible that the stars formed in their observed locations 
within the last ~ 1 Myr, but this seems less likely given 
the distribution of molecular gas in relation to the stellar 
positions and the high fraction 50%) of wTTSs in the 
association. 

The radial distribution of stars fans out to the east, and 
there are significantly more stars to the east o f the cloud 
center than the west (see Figures [2l and |3 



Moriarty- 



Schieven et al.l (|2006|) argue that LT551 isbeing eroded 
from the east by ionizing flux, perhaps from Orion, and 
the PMS stellar distribution is modest support for their 
conclusion that star formation is progressing predomi- 
nantly from east to west. 

Although we are only observing projected separations, 
the spatial distribution of embedded stars does not ap- 
pear to follow any particular prescription such as Jeans 
fragmentation. All embedded stars appear in regions of 
stellar activity sugges ting that triggering ma y be an im- 
portant process {e.g., Yokogawa et al. 2003). However, 
the fact that the next star or stellar syste m m L I 551 ap- 
pears to be forming from quiescent gas (§3.4 1 counters 
this notion. 

3.4. Future Star Formation 

A gravitationally bound, 2-3 Mq pre-protostellar core 
was recently discovered in a quiescent reg i on to the 
northwest of L1551 I RS5 ([Swift et al.||2005j [Moriarty- 
Schieven et al.[ [2006| . Follow up observations of this 



core, named Li551-MC, confirm the existence of high 
volum e density gas and reveal > O.lkms"^ infall signa- 



tures ( Swift et al.[[2006 1. The gas motions related to the 
redshiited self-absorption in the CS line profiles are as- 
sociated with build up of core mass, and it is probable 
that a star or stellar system will form in this region in 
< 1 My r. However, the existenc e of a very low luminosity 
object ( [Kauffmann et al.[[2005] embedded in L1551-MC 
has not yet been ruled out. 

Figure [?[ shows a composite image of L1551-MC made 
from all four IRAC bands. There is a deficit of emis- 
sion along the central dust lane in the 8 /im band due 
to the extinction of the Galactic background at these 
wavelengths. In the 3.6 and 4.5 fim bands there is excess 
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Fig. 7. — Spitzer IRAC composite image of L1551-MC with blue 
= 3.6 iim, green = 4.5 iim+5.8 p,Ta, and red = 8.0 fim. The position 
of peak NH3 emission from jSwift et al.| | |2005j l {white cross) and an 
ellipse enclosing the highest column density region of the map are 
shown. 



emission along the lane due to scattered light. There is 
no evidence for an embedded source within an arcminute 
of the peak NH3 emission in our IRAC data. An upper 
limit of Lboi < 0.005 L0 is set from the sensitivity of 
our 8 /im IRAC image using a bolometric temperature of 
100 K. 

Two 24 /im sources appear near L1551-MC in 
MIPS archival data (Campaign ID: 713). The 
brighter source has a flux ^ 1 Jy and is located at 
4^31™05!7,+18°13'21" (J2000). This source has no 
counterpart in any of the IRAC bands and shows a 
systematic 0"22s~"'^ motion relative to another point 
source in the MIPS field. It is likely this source 
is an asteroid. The weaker source is located at 
4i^31™07%+18°13'17" (J2000) and is coincident with an 
extended source apparent in all IRAC bands. This ob- 
ject may be either multiple sources along the line of sight 
or a background galaxy. 

The 3.6 and 4.5 /im images are used to estimate the 
total visual extinction inside and outside the central 
2' X 3'.5 region of L1551-MC shown in Figure [t] The 
mean [3.6] — [4.5] color of stars inside this oval translates 
to {Av,in) = 42™ ± 4™. The colors of background stars 
outside the oval give (Av^out) = 19™ ± 3™, consistent 
with the NIR measurements in this region. 

A low-mass star or stellar system is likely to form 
within the high column dust lane extending northwest 
from L1551 IRS5, and these new Spitzer data show that 
a low luminosity protostar has yet to form in the L1551- 
MC region. Given the limited extent of the dust lane 
(see Figure [2]) , it is unlikely that a significant number of 
stars will form in the future. However, further study is 
needed to determine the gravitational stabihty of the gas 
northwest of L1551-MC. 



4. GENERAL PROPERTIES OF THE L1551 CLOUD 
4.1. Mass and Abundances 

The visual extinction shown in Figure [2| is conve rted 
to total hydrogen column using Bohlin et al. (19781 and 
an extinction curve slope, Ry = 3.1, giving iV(Htot) = 
1.9 X lO^^Aycm^^. Summing over the entire cloud 
yields a mass of Mh ~ 120 M0. Correcting for a he- 
lium abundance, Yp = 0.25, the total mass of L1551 
^tot ~ 160 M0. This value is higher than all previous 
estimates of the ambient cloud mass made using CO iso- 
topologue emission {e.g., Snell et al. 1980 Stojimirovic 
[et al.,^2006} . CO fails to trace areas of low extinction 
due to photodissociation and also suffers saturation ef- 
fects in high column regions. These effects are expected 
to account for the discrepancies. 

Assuming local thermodynamic equilibrium (LTE) 
conditions, the optical depth of both ^^^00(1 — 0) and 
0^^0(1 - 0) can be estimate d over the central 20' x 20' 
of LI 551 dSwift et al!]|2005[ Equation 3). The maxi- 



mum "CO optical depth in L1551 is computed to be 
■''IS, max ~ 10, and the mean (ria) ^ 2. The C^^O emis- 
sion is optically thin everywhere. 

A pixel by pixel comparison of our Ay map with the 
regridded CO isotopologue maps gives a conversion be- 
tween magnitudes of visual extinction and total CO iso- 
topologue column. We assume a constant excitation tem- 
perature of 15 K for the molecul ar gas, although varia- 
tions from 9 to 25 K may exist dSwift et aL 2005[ Sto- 



jimirovic et al. 2006 



Snell 1981 r I'here is a linear 
correlation between both the r-corrected ^^CO and the 
C^^O column depth with Ay between ~ 2 and 10™. For 
Ay ?J 10™ the scatter in the relation becomes large and 
the CO derived column depths appear to saturate. Lin- 
ear fits to the column depths as a function of Ay between 
and 10™ give 

iV(i3C0) = (2.5 ± 0.1) X 10^^ (Ay - 0.8 ± 0.3) cm-^2) 
N{C^^O) = (2.6 ± 0.2) X 10" {Ay ~ 2.2 ± 0.4) cm-^3) 

with the 1 a error estimates from the fits included explic- 
itly. 

These numbers agree well with pa st studies (e.g ., 
man||1978^ 'Lada et al.||1994|). Using Bohhn et al. 

^•^CO is 1.3 ± 0.05 



Dick- 



IS 



TTT^fST^ 



the abundance ot 
between 0.8 and 10™, and the abundance of C^^O is 
1.4 ± 0.07 X 10~^ relative to the total hydrogen content. 
Using an excitation temperature of 25 K increases the 
abundance of the isotopologues by ^ 25%. 

4.2. Energy 

The gravitational energy of L1551 depends on the to- 
tal mass and the distribution of matter. A radial profile 
is computed using averages of the r-corrected ^'^CO and 
Ay maps within circular annuli centered on the C^^O 
intensity weighted mean position of L1551, 4'^31™24'^, 
+18°10'00" (J2000). Figure [8] shows this column den- 
sity profile from ^ 0.02 to 2pc. The error bars for each 
value are determined from the rms deviation of pixel val- 
ues and the number of pixels in each annulus. The errors 
derived for the ^^CO data lie within the symbols. The 
profile is fiattened in the inner region of L1551 and be- 
comes progressively steeper out to the edge of the cloud, 
^^cloud = 0.9pc. 
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Fig. 8. — Normalized column density profile of the L1551 dark 
cloud derived from both r-corrected ^"^CO emission (open circles) 
and near infrared extinction (filled circles). The maximum aver- 
aged column density is (A'"(Htot))max = 2.2 X 10^^ cm~^. 



For a density profile 



Po 

Po (Rcorc/rY 



< r < i?core 
-^corc — ^ — -^cloud 



with Rcorc = 0.1 pc, i^cioud = 0.9 pc, 
160 M0, the gravitational energy 



rpcoTC 
grav 



3 GM, 
5 



2 

core 



Rc 



1 



loud 



and Mu 



^core) 



Rc 



(4) 



(5) 



is — 8.7x 10 ergs, where po = 1.9 x 10^ gcm^'' is found 
by normalizing the total mass to 160 M©, and Mcore = 
4/3 Trpo^core- Upper limit to i?grav is given by the 
gravitational energy of a truncated, singular isothermal 
sphere (SIS) with the total mass and radius of L1551, 
EgSv = ~GAP/R = -2.5 X lO-^^ergs. In Hght of these 
calculations, we use a gravitational energy of -Egrav = 
— 1 X 10"*^ ergs in further analyses with a factor of < 2 
uncertainty. 

The turbulent velocity width derived from the compos- 
ite C^*0 spectrum (see §[7| is utmb = 0.52 kms^^. The 
one-dimensional thermal width of the molecular gas 



C^tho 



(6) 



is 0.23kms-i for T = 15 K and to = 2.3 toh- Adding 
in quadrature the turbulent and thermal contributions 
gives the three-dimensional velocity dispersion of the gas 
in L1551, (Tv — 0.65 kms^^. The total kinetic energy is 
thus 6.8 X 10*** ergs, very near what the virial theorem 
predicts given our estimation of -Egrav • 

4.3. Timescales and the Jeans Criterion 

Approximating L1551 as a homogeneous sphere with 
density (n) « 1000 cm~'^ and mean particle mass m = 

2.3 Toh, the free-fall timescale tg = [37r/(32 Gp)]^^^ — 
1.1 Myr. The dynamical timescale of the cloud is taken 
to be equal to the free-fall timescale, idyn = iff- The 
sound crossing time in L1551 is estimated as icross = 
2i?cioud/cs = 7.7 Myr for = {kT/rfif'^ and Tk = 15 K. 
This means that L1551 cannot be supported by thermal 
pressure alone. A varying temperature structure could 



alter this number significantly. However, a kinetic tem- 
perature of > 700 K would be needed for icross = iff- 

Another assessment of the balance between t hermal 
press ure and gravity is given by the Jeans criterion ( [Jeans 
1961 1. For a uniform, isothermal gas of density p, pertur- 

bations exceeding a length scale of Aj = {tvc^ /Gp) will 
be unstable to gravitational collapse. This length scale in 
L1551 is Aj « 0.8 pc — about half the cloud diameter — 
meaning it is Jeans unstable. Therefore other mecha- 
nisms, such as magnetic fields or turbulence, must con- 
tribute to the overall gravitational stability of the cloud 
to prevent a global collapse. 

Unfortunatelj/, little is known about the magnetic field 
in L1551 (see §p|. However in the next few sections, a 
closer look at tne molecular gas in L1551 gives insight 
into sources and dissipation of kinetic energy that are 
keeping this cloud close to virial balance. 

5. THE MOLECULAR GAS AT HIGH-RESOLUTION 

The ubiquity of outflow from young stellar objects and 
the presence of ^ 35 proto- and pre-main sequence stars 
associated with LI 551 imply a complex relationship be- 
tween the stars and the gas in L1551 over its lifetime. 
The interplay between young stars and the gas from 
which they formed is best revealed in our fully sampled, 
high-resolution images of ^^CO and C^^O emission. 

5.1. Velocity Integrated Emission 
Figures [9] and 10 show the velocity integrated ^"^CO 



and C^*0 emission from L1551, respectively. These 
maps have been created by summing ^T^AV^h over 
the line profile, where AVch = 0.13 kms~^ is the channel 
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Fig. 9. — Velocity integrated ^'^CO(l — 0) emission of the com- 
bined interferometric and single-dish data. The image maximum is 
20.6Kkms~^. The axes are labeled in reference to the position of 
L1551 IRS5, 4'>31™34!1, -|-18°08'04" (J2000). Sources are shown 
as crosses and the locations of selected HH objects are marked with 
triangles. The thin lines represent the jet directions from the em- 
bedded sources and the resolution of the inner and outer regions 
of the map are shown at the bottom left. 
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Fig. 10. — Same as Figure |9| but for the velocity integrated 
C^*0(1 — 0) emission. The maximum pixel value in the map is 
5.5Kkms~^ at the position of L1551 IRS5, but the linear stretch 
has been adjusted to bring out lower level features. Note the sky 
coverage for the C^*0 data is not exactly the same as the ^^CO 
data. 



width. The high-resokition data are shown within the 
dashed contour, while outside this boundary there are 
only single-dish data. The C^^O data cover a smaller 
range in right ascension but a slightly larger range in 
declination. 

The ^"^CO map shows a smooth overall distribution of 
emission with clearly delineated features around regions 
of known stellar energetics. Two caliper-like arcs of emis- 
sion to the southwest of L1551 IRS5 trace the boundary 
of the southwest cavity from the well known LI 551 IRS 5 
bipolar outflow (jSnell et al.||1980| [Draper et al.||1985 



Moriarty- Schieven et al.|1987 ! Moria rty-Schieven &: Snel 



1988 Moriarty-Schieven et al.|,2006i |Stojimirovic et al. 
200?! . Northeast of L1551 NE, a sharp edge in the ^•='00 



map creates a par abolic shape that is symmetric around 
the L1551 NE jet ( [Reipurth et al.|2000[ ), though the lower 
arm is less pronounced than the upper. There is a deficit 
of ^^CO emission within this feature that correlates with 
a low extinction derived from background stars (see Fig- 
ure[2]), and we refer to this region henceforth as the north- 
east cavity. The gas around XZ / HL Tau shows clear si gns 



of disruption, likely due to jets (Mundt et al. 1990 ) as 



well as less coUi mated outflow seen in the region (Krist 



|eral.||1999| JCoffey et al.||2004[ [Welch et al.|[2000| . 

The O emission is more clumpy than the ^^CO 
emission and follows more closely the distribution of Ay . 
Most features seen in C-'^^O correlate with features seen 
in -'^'^CO, but there are a few notable differences. The 
northwest region of the -'^'^CO map near HH265 is rel- 
atively featureless, while the smooth lane of extinction 
that encompasses L1551-MC is noticeable in C^^O. A lin- 



ear feature seen in C^^O extending southwest from the 
XZ/HL Tau region has no counterpart in the ^"^CO map, 
and the southern edge of the southwest outflow lobe from 
L1551 IRS5 is not detected in C^^O. 

5.2. ^'^CO Velocity Centroids and Widths 

Figure [TT] shows the first and second "moment" maps 
of the ^'^CO emission. Gaussian fitting of the ^^CO line 
profile in each independent, 10" x 10" pixel produced a 
velocity centroid, Veen (left panel) and velocity full width 
at half maximum, AV = 2^/2 In 2a (right panel) used to 
construct these maps. The Gaussian shape is a decent 
representation of the ^^CO line profile with ^ 85% of the 
fits having values less than 2. 

The intensity weighted mean velocity of the ^"^CO emis- 
sion is 6.7kms~^. Some of the smallest values of V^ccn he 
along the border of the southwest lobe of the L1551 IRS5 
fiow. The diffuse emission in the northeast cavity also 
has v ery blue central velocities likely due to outflow (see 



5.3 1. The red velocities at the w estern edge of the map 



are related to the L1551W flow (Pound fc Bally 1991 



Moriarty-Schieven fc Wannier[[l991[ [Stojimirovic et al 
2006 1. Variations at the 0.2 kms~^ level are found across 
the inap that loosely correlate with features in Figure [9] 
The regions where young stars are injecting energy into 
the cloud can be cle arly identified in the velocity width 
(&). The peak hue widths exceed 
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map shown in Figure 

2kms~^ in a very thm (< 1600 AU in some sections) 
rim around the southwest lobe of L1551 IRS5 and in the 
northeast cavity. High line widths are also seen around 
HL/XZTau and in the L1551W fiow region. The mean 
^•^CO line width in L1551 is 0.96 km s^^ with a variance 
of ~ 0.05kms-i. 
The ^■^CO velocity widths in Figure 



11 



b) are a proxy 

for the kinetic energy in the gas along each line of sight. 
The velocity widths are only significantly above average 
in the regions directly associated with the energetics from 
young sources. Small spatial regions containing high line 
widths like the border of the southwest lobe have small 
diffusion timescales and suggest that these structures are 
the result of current outflow activity. Stellar feedback is 
the most significant source of energy input in the cloud, 
and the only source noticeable in our -'^^CO data. 

5.3. The Spatial and Kinematic Distribution of ^■^CO 

The nature and relevance of the features seen in Fig- 
ures [9] and [n] are further revealed in the ^^CO chan- 
nel maps. Figure 
velocity intervals c 
tures in the 



12 shows the structure of L1551 at 



;hosen to highlight the numerous fea- 
^■^CO data cube. The central velocity of 
each map is increasingly blueward of the line center in 
images (a)-(d) and increasingly redward in images (e)- 
(/i). Each channel map is depicted with a linear trans- 
fer function for the gray scale shown to the right of the 
plot adjusted to cover all integrated intensities from the 
2.5a level to the maximum value in each map. These 
values are; (a) 0.40-7.4 K km s"!; (&) 0.33-7.6Kkms-i; 
(c) 0.35-4.1 K km s-i; (rf) 0.30-3.4K kms'^; (e) 
0.45-6. 5Kkms-i; (/) 0.45-3.2Kkms-i; (5) 0.33- 
3.3Kkms-i; (h) 0.35-2.1 Kkms"!. 

5.3.1. Synopsis of Channel Maps 

In the line core channel maps (a) and (e) much of 
the emission is optically thick, and the disruption of the 
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Fig. 11. — The first and second "moment" maps of the ^^CO emission in L1551 derived from Gaussian fitting. Sources and HH objects 
from Figurejgjare shown as well as the jet directions from the embedded sources. Spectra with '^^CO line peaks below 3 a produce unreliable 
fits and are masked in white. The dashed black line encloses the high-resolution data. 



ambient gas by outflows creates the most apparent fea- 
tures. In channel map (a), the northeast cavity appears 
as an elliptical hole with its major axis aligned with the 
LI 551 NE jet, and the southwest cavity is outlined by 
a thin parabolic arc. The limb brightened shell arou nd 
XZTau investigated in depth by Welch et al. (20001, is 
apparent along with a few linear structures to the east 
of the XZ/HL Tau region. 

In channel map (&), the diffuse, blueshifted emission in 
the northeast cavity responsible for the small T4en values 
in this region [Figure 111] (a)] is clearly visible. The bright 
limb around the southwest lobe becomes a dominant fea- 
ture at these velocities, and the brightest emission em- 
anates from a thin wall of compressed molecular gas just 
downstream from HH 102. Low level emission present 
near the "mouth" of the southwest cavity lies along the 
axis of the L1551 NE jet. 

Moving further into the blueshifted line wing in chan- 
nel maps (c) and (d), the bright rim of the south- 
west cavity persists. The emission within the main 
L1551 IRS5 flow moves downstream but remains com- 
pact and roughly aligned with the L1551 NE jet. Knots 
of emission appear around the XZ/HL Tau region. 

On the red side of the ^^CO line core shown in channel 
map ( e) , three regions of proto-stellar disruption are visi- 
ble, the southwest cavity associated with the L1551 IRS5 
outflow, the northeast cavity, and the XZ/HL Tau region. 
The southwest cavity is outlined by a complex pattern, 
while the northeast cavity shows a sharp upper bound- 
ary. The lower edge to the XZ Tau bubble is prominent 
while the emission to the north of XZ Tau appears frag- 
mented. 

In channel map (J) the emission is broken up into a 



network of thin filaments. The lower rim of the XZ Tau 
bubble is clearly visible and extends to the position of 
HH30IRS. The upper limb of the northeast cavity has a 
filamentary appearance at these velocities. Northeast of 
L1551 IRS5, between L1551 IRS5 and L1551 NE, there 
are a series of arcs and knots of emission. Perhaps the 
most interesting is the arc that intersects the position of 
L1551 IRS5 and may trace the boundary of its red out- 
fiow lobe. The bright emission west of HH 10 2 remains 
a prominent feature, and the L1551W flow ([Pound k. 
Bally|199T Moriarty-Schieven fc Wannier|1991"| appears 
filamentary in the high resolution part of the map. 

Continuing further into the red-shifted velocities in 
channel maps (g) and (/) the arc intersecting L1551 IRS5 
persists but moves gradually northeast and fades with in- 
creasing velocity. A few knots of emission remain in the 
highest redshifted velocity map likely associated with the 
red lobe of the main L1551 IRS5 flow. 

5.3.2. Discussion 

The high-resolution ^"^CO data show a wealth of struc- 
ture never before seen in this well-studied region. It is 
clear from this map that energy from embedded sources 
has accelerated ambient cloud gas. Some of this gas re- 
sides in the line wing emission. A vast majority of the 
features in the optically thin ^'^CO line wings appear as 
thin arcs of emission. A lo w filling fac t or ha s been de- 
duced for the L1551 flow by Snell et al. (19841 from low- 
resolution data, and the structure seen in the "'^'^CO line 
wings of Figure 12 is consistent with their value of '--^ 0.3. 



The diffusion timescale of distinct features seen in the 
channel maps can be estimated by dividing the size of the 
feature by the width in velocity space, tdiff ~ Aa;/Au. 
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Fig. 12. — Channel maps of ^''CO emission in L1551. The positions of all known embedded sources (crosses), selected HH objects 
(triangles), and jet directions from sources (thin lines) are shown. See the text for a detailed description of these maps. 
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Typically, Aa; ~ 0.05 pc (see §7.2) and Av ^ Ikms^^, 
giving tdis ~ 10^ years. The coherence in both position 
and velocity space of the filamentary structures seen in 
the ""^^CO line wings along with the short implied diffu- 
sion timescales leads us to interpret them as the shocked 
gas boundaries separating quiescen t gas from outflow mo - 



tions from young sources (see also Barsony et al.||1993 l 
The velocity widths of these features do not extend much 
beyond the escape velocity of the cloud, ^ 2kms~"'^, and 
will likely continue to decelerate into the ambient gas. 
The large Reynolds numbers of these flows imply that 
they are highly turbulent. 

6. STELLAR FEEDBACK 

The high-resolution maps of §[5] show that energy is 
being injected into the ambient cloud from where known 
embedded sources exist. We look at three particular re- 
gions that show clear signs of stellar feedback — the south- 
west lobe of the L1551 IRS5 outflow, the northeast cavity, 
and the XZ /HL Tau region — to gain insight into the how 
stellar feedback is affecting the molecular environment of 
L1551. 

6.1. The Expanding Shell from the L1551 IRS5 
Outflow 

The blueshifted, or southwest, lobe of the L1551 IRS5 
outflow a p pe ars a s a thin, limb-brightened shell in Fig- 
The position- velocity {p-v) cuts shown 



and 12 



ures[9j[10 _^ 

in Figure 13 oriented perpendicular to the L1551 IRS5 
outflow axis reveal the kinematic nature of this feature. 

Figure 14 shows the emission in p-v cuts (a)-(rf) in 
order of increasing distance from L1551 IRS5. Along cut 
(a) a full arc of emission blueward of the line core can 
be seen. This p-v signature is what is expected from an 
expanding shell with a velocity of ~ 2kms~^. Moving 
further away from LI 551 IRS5 in cuts (&) and (c), the 
shell walls appear spatially thin, ~ 0.02 pc, and have 
large velocity widths. Beyond velocities 3-4kms~^ from 
the line core the shell boundaries become discontinuous 
in p-v space, and in cut (rf) no significant amount of 
emission appears at high velocity. This suggests that the 
expanding shell has burst out the front side of L1551, 
in accord with previous interpr etations of observational 



data ( e.g., Draper et al. | 1 985; Moriarty-Schieven & Snell 
2000D . The dynamical age of this 
4'. 5, and expansion velocity 
years 



1988| [White et al| 
feature from its diameter , 
^ 2kms~^, is roughly 10^' 



The symmetry of this expanding shell is not aligned 
with then kno wn direction of the L1551 IRS5 jet (Mundt 
fc Fried|[T983 l as seen in Figure |l3l Yet the shell wall is 
spatially thin, ~ 0.02 pc, and has a large spread in ve- 
locity, ^ 3kms~^, implying a short dispersion timcscale, 
Ax/Av ~ 10'' years. High velocity H i seen in t he outflow 



lobes of L1551 IRS5 (Giovanardi et al. 2000[ ) may con- 
tribute to the acceleration of this shell and its spatially 
thin structure. 

Summing the emission in the blue lobe at LSR veloci- 
ties less than 6 kms~^ gives a total mass of ~ 3 Mq. This 
is taken to be a lower limit to the total shell mass since 
the shell also contains a component at the line-of-sight 
velocity of the ambient gas (see Figure 14 1 . The energy in 
the expanding shell estimated by converting ^^CO emis- 

T4cnP in each 
< 10^'^ ergs. 



c 



^5 
<3 




-2 -4 
Aq (arcmin) 

Fig. 13. — Velocity integrated ^''CO emission as in Figuretolwith 
the position-velocity cuts through the blue lobe of the L155iIRS5 
outflow overlaid (a-d). The positions of L1551 IRS5 and promi- 
nent HH objects (triangles) are labeled, and the direction of the 
L1551 IRS5 jet is shown (arrow). 
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velocity channel less than 6 km s ^ is 



Fig. 14.— i^co emission in the position-velocity cuts across 
the blue lobe of the L1551 IRS5 outflow at a position angle of 
135°, roughly perpendicular to the L1551 IRS5 outflow axis. The 
projected distances of the cuts from L1551 IRS5 are ^ 0' 5. 1 '. l'.5, 
and 2' for (a), (b), (c), and (d), respectively (see Figure [l3| . The 
gray scale follows a linear transfer function from 0.3-7 R and the 
contour levels are 0.9, 1.5, and 2.1 K. The zero position roughly 
delineates the outflow axis of symmetry. 



The total mass traced by molecular emission in the en- 
tire blue lobe of the L1551 IR S5 outflow is higher than 
the value derived for the shell (j Stojimirovic et ar]|2006 l. 
However, the high velocity gas concentrated toward the 
center of the outflow lobe seen in CO will most likely 
end up in the diffuse interstellar medium rather than 
the L1551 cloud. The existence of HH objects beyond 
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the edge of the cloud (see Devine et al. 1999 Moriarty- 
Schieven et al.|2006 Stojimirovic et al.|2006p and the re- 
tlection nebulosity emanatine from L1551 IK 



bulosity emanating from L1551 iiiSS (Draper 

et al.|[l985 l is further evidence that a fraction of the ac- 
celerated ambient gas is being lost from the cloud. 

An estimate for the mass of excavated material can 
be made by calculating the decrement of ^^CO emis- 
sion inside the southwest cavity relative to the ambient 
cloud. Using Equation [2] the total mass excavated from 
the southwest lobe is of order 2 Mq . 

6.2. The Northeast Cavity 

The lack of CO isotopologue emission and extinction 
of background starlight in the region to the northeast 
of L1551 NE, the shar p b oundary around this deficit 
seen in Fig ures M and 
eTaLlpOOe] Figure 11) 



12 (also see Moriarty-Schieven 
_ and the pres ence of HH 286 be- 
yond the cloud edge in this d irection (De vine et al.|1999 
Moriaxty-Schieven et al. 2006 ; Stojimirovic et al. |2006[ 
are all evidence that mass has been excavated from this 
region by stellar energetics. This feature has been seen 
in previous ly published data (see Moriarty-Schieven & 
Snell|1988| Figure 4; Pound & Bally 1991, Figure 4; Sto- 
jimirovic et al. 2006, Figure 3), but has never been in- 
terpreted. 



In Figure 15 our deep NIR imaging is shown as a JHK 
composite image overlaid with ^^CO contours. NIR neb- 
ulosity on the southwest side of LI 551 NE is conical with 
an opening angle of ^ 60° and i s symmetric wit h the 
jet ax is of - 243° ( |Devine et al.|[T999l |Rcipurth et al 



2000 



On the northeast side, the nebulosity follows the 
parabolic shape of the ^■^CO emission and is brighter 
along the northern arm. This nebulosity extends from 
the LI 551 NE source position and has a geometry that 
implies alignment with the jet axis. 
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Fig. 15. — Deep JHK composite image of the L1551 NE re- 
gion overlaid with contours of ^^CO emission integrated from 6.5- 
7.9kms~^ (LSR). The linear stretch on all NIR bands is from —a 
to 20(T, where cr is the rms of the image, and the contour levels 
are (n -|- 2) Kkms~^, where n is an integer. L1551 NE is the red 
source to the lower right in the image and the position of HH 262 
is labeled {white triangle). 



IIII262 lies within the northeast cavity and has com- 
ponents with a wide range of radial velocities and proper 
motions suggesting that the jets from both L1551 IRS5 
and L1551 NE may be interacting w ith this region 
(Devi ne et ar]|1999| [Topez et al.|ll998[ ). Faint K-hand 
emission seen ^ 1'.5 to the northwest of the nominal 
IIII262 position in Figur e [TS] is likely due to shocked 
II2 associated with GII9 ( Graham fc Heyer||1990 Lopez 
ei al. 1998). HH286 lies further out of the L1551 cloud 
and atsb origin ates from either L1551 IRS5 or L1551 NE 



pevine et al |1999^ . The red lobe of the L1551 IRS5 
molecular outfiow overla ps with the northeast cavity (see 
Stojimirovic et al. 2006 Figure 7), but the axis of that 



tiow lies along the northern arm of the cavity making 
L1551 IRS5 unlikely responsible for excavating the mass 
from the northeast cavity. 

We therefore attribute the northeast cavity feature 
to L1551 NE while, similar to the blue lobe of the 
L1551 IRS5 outflow, it contains features from both 

th 



L1551 IRS5 and L1551 NE dDevine et al. 



dynamical time for the northeast cavity 
projected size and a velocity spread of 



1999]). The 

ased on its 



3kms is 



10 



4-5 



years. The mass of gas excavated from the re- 
gion is estimated to be 2-3 Mq from the deficit of ^^CO 
and Av in the cavity relative to the average column 
depth of the cloud around the cavity. 

6.3. The XZ/HL Tau Region 



Th e XZ/HL Tau region is highly active (Mundt et al 
1990[ ). XZTau is known to have poorly collimated out- 
bursts on several year timescales (Krist et al. 1999), and 
a limb-brig htened shell of ^ ^CO is centered roughly on its 
position ([Welch et al.|2000[ ). In Figures [T2|; a) and (e) the 
acceleration of the anibient gas from around XZ /HL Tau 
is clearly evident. There is also a circular region devoid 
of ^■^CO emission in channel map (a) centered on the 
position of LkHa 358 that may be evidence for a weak 
wind. 



The expanding half shell seen by Welch et al. (2000) is 
evident in channel maps (a) and (e) and extends t o th e 
position of HH30IRS in map (/). Also in Figure p];/) 
there are two arcs curving around the east side of XZTau 
that are disjoint from the lower rim of the XZ Tau shell. 
The near one is smooth and filamentary, and the far 
one is clumpy. These features persist for multiple chan- 
nels and move in unison away from the XZ/HL Tau re- 
gion with decreasing velocity. This is consistent with 
gas expanding into the ambient cloud as it is accelerated 
from the XZ/HL Tau region. The expansion velocity is 
^ 0.7kms~^ and the expansion age is < 10^ years as- 
suming this scenario. A simple estimate for the time 
between the ejection events is ~ lO'* years. It is possible 
that these are molecular signatures o f episodic outbursts 



similar to those s een in visible light (Coffey et al. 2004 



On the 

ures 



d1[l999l. 



Krist et 

ueshifted side of the ^'^CO line core in Fig- 
12 'a) and (6) a peculiar linear feature located ^ 2' 



east oiXZ Tau extends in a northeast-southwest direc- 
tion . This feature persists over ~ 2kms~^ and does not 



vary in position over those velocities. Figure 16 shows a 
composite JHK image of the XZ /HL Tau region overlaid 
with ^-^CO contours. This linear feature seen in ^-^CO 
coincides with NIR nebulosity of the same morphology 
that is likely due to reflected light from XZ/HL Tau. 
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Fig. 16. — Same as Figure [Ts] but for the XZ/HLTau region. 
Tlie ^^CO map has been integrated over velocities between 5.4 and 
7.7kms~^ and the contour levels are (5 + 1.5n) Kkms~^. 



This hncar feature is also visible in Ha emission fDevine 
leF al. 1999J , and likely delineates the rim of the cav- 
ity vacated by the energetics from XZ Tau, HL Tau, and 
HH30IRS (see |Mundt etal]|l990l ) . 

The ^'^CO emission shown m li'igure 16 also traces well 
other diffuse NIR nebulosity. The nebulosity immedi- 
ately below XZ Tau and the nebulosity stretching north- 
ward from HL Tau also may be delineating a boundary 
between high and low extinction regions in L1551. Un- 
fortunately, not enough background stars are seen in our 
deep NIR images to do an extinction analysis in this 
region. Howev er, these data are c onsistent with the in- 



terpretation of Welch et al. ( 2000 1 that the HL Tau lies 
at the boundary of a region excavated by XZ Tau. 

Stars in this region of L1551 must have destroyed part 
of the cloud both to expose XZ Tau along our line of 
sight and to create the features seen in ^■^CO and the 
NIR. However, the modest diminution of ^"^CO column 
depth seen in our figures translate to only ~ 0.1 Mq of 
excavated matter from the XZ/HLTau region. 

6.4. Summary of Stellar Feedback 

A number of new features are revealed in our high- 
resolution maps that show clear signs of interaction be- 
tween outflow from young stars and the ambient cloud 
from which they formed. The ^'^CO line wing emission 
traces the interaction between gas accelerated by embed- 
ded sources and the ambient cloud. At our sensitivity 
limit, this emission does not extend to velocities much 
beyond the escape velocity of the cloud and is therefore 
a good tracer of the energy that will remain in the cloud 
and ultimately contribute to its turbulent velocity field. 

Using the closely Gaussian C^^O composite spectrum 
to characterize the quiescent gas, the ^'^CO line wings 
are defined as V^ed > 7.5kms~^, and Vbiuo < 5.9kms~^ 
(LSR). The energy in the red and blue wings are '-^ 1.5 
and 2.0 x lO^'^ergs, respectively. If isotropy of outfiow 
motion is assumed, the total energy in the ^■^CO line 



wings becomes ~ 10**^ ergs. 

This is about a factor of 5 lower than recent estin iates 
for the total outflow energy in L1551 using ^^CO (|Sto- 



jimirovic et al. 2006) which traces both the bound and 
unbound outflow gas. All major outflow regions show 
clear evidence that gas is being accelerated beyond the 
gravitational confines of the cloud, and at least 4-5 Mq 
of mass has been excavated from L1551 by the present 
day embedded stars. This suggests that L1551 was more 
massive at the inception of star formation by perhaps 
60%. 

7. TURBULENCE 

While our ^'^CO data show a wealth of structure in the 
line wings due to the energetics from young embedded 
sources, the C^^O emission has a closely Gaussian profile 
at our sensitivity level and is optically thin throughout 
the map. The C^^O emission is thus better suited to 
probe the turbulent velocity and density fields in L1551. 

The one-dimensional velocity dispersion of the compos- 
ite C^^O spectrum is 0.31 kms^^. Subtracting in quadra- 
ture the thermal contribution to the measured width us- 
ing Equation |6] with m = 30 mn and then multiplying 
by v^, the three-dimensional turbulent velocity disper- 
sion (Tturb ~ 0.52 km s~^. The total turbulent energy in 
L1551 is estimated to be l/2Mtotcr^^^b = 4.3 x lO'^'^ergs, 
and £:turb ~ 0.4£;gi.av 

7.1. Spectrum of Fluctuations 

A complete statistical description of fully developed 
turbulence is obtained in the scaling relations of the 
density and velocity fields. Under ideal circumstances, 
energy is conserved as it cascades from large to small 
scales at a rate proportional to v^{L)/L, implying a self- 
similar velocity dis persion th at scales as v{L) oc , with 
7=1/3 (Kolmogorov 19411. This line width-size rela- 
tion translates to an energy spectrum E{k) oc , with 
P — 5/3, and a power spectrum of Pv{k) (x k"", with 
n = —{j3 + D—l) and D being the number of dimen- 
sions that the power spectrum is computed ove r. Other 
turbulence models, suc h as shock dominate d (Burgers 
1974 1 and multi- fractal ( Boldyrev et al.|2002 1 , also show 
distinct scaling relations. 

Unfortunately, observers can only measure two dimen- 
sional distributions of velocity and density on the plane 
of the sky, and many different techniques are used to de- 
duce the three-dimensional distributions from the avail- 
able information. Line widths and velocity centroids 
are a direct measure of the velocity field and have been 
used widely to in fer the charac t eristics of the underlying 



turbulence {e.g., (Larson||1981| iScalo 1984 Dickman fc 



Kleiner 1985: Kleiner & DickmaK p^85| [Miesch Ijally 
|1994j . Other meth ods including the spectr al correla- 
" ' ' . - , prin ciple 

2002|, A- 



tion function (SCF; Rosolowsky et al 



compone nt analysis (PC A; Brunt & 



I 



1999 



eyer 



variance ( Stutzki et al.|1998t, and veloci ty "cha nnel anal 



ysis (VGA; Lazarian & Pogosyan 2000 , 2004 1 have had 
varying degrees of success. 

7.1.1. Velocity Channel Analysis 

Density and velocity fluctuations are not independent 
within the framework of homogeneous, isotropic turbu- 
lence and can theoretically be discerned through observa- 
tion. The velocity channel analysis (VGA) technique has 
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Fig. 17. — The annular average of the two dimensional density 
power spectrum of C^^O emission as a function of wave number. 
The corresponding scale in the map is plotted on the upper axis. 
The best fit power-law slope for Pp^ik) is —2.8 and is shown as 
the line traversing the data. 



been develo ped toward this end ([Lazarian & Pogosyan 



2000 2004) and has been used to successfully separate 



the velocity and density contributions to th e power spec- 



trum of H I in the Small Magellanic Cloud (jStanimirovic 
fc Lazarian||2001| . 



i'he intensity fluctuations within a spectral line map 
integrated over the full extent of the line profile are dom- 
inated by the density fluctuations such that the index of 
the two-dimensional power spectrum of intensity fluctu- 
ations equals the index of the three-dimensional density 
power spectrum 



P, 



2D, thick 



(fc) cx Pl°{k) cx k 



(7) 



Intensity fluctuations within a thin slice in velocity 
space have contributions from both the velocity and 
density fluctuations in proportions that depends on the 
steepness of the three-dimensional density spectrum 



52D,thii 



(fc) cx 



n > 
n < 



(8) 



Here a is the power-law ind ex of t he second-order veloc- 
ity structure function (see §7.1.31. The formal criterion 
for a slice to be thin is that the dispersion of turbulent 
velocities on the scales studied be larger than the veloc- 
ity slice thickness. The thin slice regime is not attainable 
for sub-sonic turbulence. 



Figure 17 shows the power spectrum of the velocity 
integrated C^^O within the high-resolution region out- 
lined with the dashed line in Figure [TO] The velocity 
integrated map was flrst buffered on all sides with zeros 
to minimize aliasing effects and was then Fourier trans- 
formed using a Fast Fourier Transform (FFT) algorithm 
and squared. Circular annuli are averaged to obtain the 
mean value of the power spectrum in each annulus, and 
the dispersion of values within each annulus determine 
the error bars. 



The power spectrum of Figure 17 follows a power-law 
form over the full range of scales sampled. A non-linear 
least squares flt to the power spectrum returns a slope 
of n = —2.8 ±0.1, characterizing the density fluctuations 
as "shallow," n > —3. The error here is derived from the 
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Fig. 18. — The derived power-law index for Pp^{k), the two- 
dimensional power spectrum of density fluctuations, as a function 
of the velocity width of the integrated map. The decrease in in- 
dex as a function of increasing velocity width is expected, and the 
flattening out at large Ai> signifies the density dominated regime. 

variance of power spectrum indices flt to subsets of the 
C^^O data. 

Figure[T8]shows the slope of the two-dimensional power 
spectrum as a function of velocity width. For each width. 
Aw, the C^^O map was integrated in velocity across the 
necessary number of channels centered around the cen- 
tral channel in the line proflle, Vlsr = 6.7kms~^. A 
power spectrum was then obtained in the same manner 
as described above. Shifting the window of channels by 
a single channel width in both directions gives two more 
measurements of the slope. The mean of the three values 
for n2D is represented as a flUed circle in Figure [18] and 
the spread in these three values determines the error bar. 

As the channel width increases, the slope of the two- 
dimensional power spectrum monotonically decreases. 
This is the behavior predicted by the VCA analysis. 
The trend levels out for channel widths greater than 
^ 0.6kms~^ and likely signifles the transition into the 
density dominated regime. There is no sign of this trend 
tapering off for small Aw and it is possible that this 
trend will continue to the thermal width of the C^'^O 
line, 0.06 kms^^ (Lazarian, A. 2006, private communica- 
tion). 

Assuming that the value of n2D = —2.44 is repre- 
sentative of the slope of Pp^{k) in the thin regime, a 
value for the velocity spectrum index can be attained 
using Equation |8] a = 0.72 ±0.1. This implies a 
three-dimensional power spectrum of velocity fluctua- 
tions, P^^{k) cx fc~'^-^^, near the value predicted by the 
Kolmogorov formalism, —11/3. 

7.1.2. Line Width-Size Relationship 

The VCA results imply an energy spectrum index, 
P — 1.72 and a line width-size index of 7 = 0.36. The line 
width size relationship for our C^^O data are shown in 
Figure [19] The one-dimensional line widths in Figure [19] 
are computed by taking an intensity weighted average 
of all independent pixels in the C^^O map downgraded 
in resolution to scale L. The average line width is then 
corrected for cloud depth effects by subtracting the dis- 
persion in line centroids at the given map resolution (see 
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Fig. 19. — The line width-size relation of the C'^^O gas in L1551. 
The sound speed for molecular gas with solar helium abundance at 
a kinetic temperature of 15 K is shown with the arrow. 
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The derived trend in line 
width is characterized by 7 = 0.07 ± 0.01. 

This derivation of the line width-size relationship is 
fundamentally different than the standard method used 
where regions are identified as individual clouds or cores 
based on prescribed criteria and their sizes and line 
widths are then determined for each separate region (e.g . , 
Larson|198l| [Solomon et al.|1987| [Fuller fc Myers||1992[ ). 



i'he observed power-law form of the line width-size re 
lationship is general ly thought to arise from interstel 



lar turbulence {e.g., Larson 19811, but the relationship 



formed in this manner does not contain any direct infor 
mation regarding the nature of turbulence within indi 
vidual clouds. 

In L1551, the line width-size relationship is very shal- 
low. This result is inconsistent with the results from the 
VGA analysis assuming that the turbulence is homoge- 
neous and isotropic. However, the derived line width-size 
relationship implies straightforwardly that the observed 
line widths are originating on small spatial scales which 
may indicate the presence of ongoing feedback (Naka- 
mura et al.||2006[ ). 



7.1.3. The Structure Function 

The velocity structure function measures velocity dif- 
ferences as a function of projected separation and is a 
common statistic used to characterize turbulent motions 



a power-law form, S'v(t) oc |t|". The index, a is related 
to the index of the energy spectrum and line width-size 
relation, a = /? — 1 and a = 2j. 

The form of the second-order velocity structure func- 
tion is 

(9) 



5v(r) = 5](K(r)-«c(r + x)]^ 



where the position, r, and lag (i.e., projected separa- 
tion), T, are two-dimensional functions on the plane of 
the sky, and the angled brackets signify an average over 



all pairs. The velocity centroid is 

j:T{r,v)vir) 



Vc{r) 



(10) 



where T is the antenna temperature and the sum is over 
velocity. The structure function of Equation [9] normal- 
ized by the variance of velocity centroids across the map. 



N 



(11) 



is denoted S* (t) , where N is the number of pixels in the 
sum and 

(12) 

is the mean velocity centroid of the map. The normalized 
structure function has an algebraic relationship to the 
spatial autocorrelation function and approach es values 
of 2 for uncorrelated spectra (see Miesch & Bally 1994). 



Detailed numerical tests have shown that the velocity 
structure function can be successfully used to retrieve 
the scaling properties of a three-dimensional velocity field 
from a data cube of position-position- velocity for turbu- 
lence with sonic Mach number Ms < 2.5 and in regions 
that are not d ominated by density fluctuations (Brunt] 
fc Heyer|[2002[ [Esquivel fc Lazarian||2005[ ). A straight- 
forward way to compare the relative importance of the 
density fluctuations is through the second order density 
structure function. 



where 



S,{T) = alY.{[Iir)-Iir + T)Y 
Iir)^J2^{r,v). 



(13) 



(14) 



The sum in Equation [14] is over velocity, and CTv is the 
velocity dispersion of the gas. To directly compare Equa- 
tion ]9] and Equation |13[ "unnormalized" velocity cen- 
troids. 



T{r,v)v{r) 



(15) 
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Fig. 20. — Normalized second-order velocity structure function 
of the C^^O emi ssion in L1551 within the high-resolution region 
outlined in Figure [To] The gray shaded region indicates the allowed 
values of S*{t) within the computed error estimates. The best 
power-law fit for r < 0.2 pc is shown. 
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replace the velocity centroids of Equation |9] The struc- 
ture function of the unnormalized velocity centroids is 
denoted Sv{t), and it has been determined by numeri- 
cal simulation that the velocity structure function accu- 
rately traces the turbulent velocity statistics for Sv{t) » 
Spir) jEsquivel fc Lazarian|[2005 1 . 

A value tor the velocity centroid was computed at each 
independent pixel according to Equation |10| Only chan- 
nels with signal greater than 2.5 a within the full width 
between first nulls of the composite spectrum were cho- 
sen as part of the sum, and only spectra with 3 or 
more channels meeting this criterion were considered, 
the rest were masked. The normalized structure func- 
tion, S*{t), is then computed using /i — 6.69 k ms~^ and 
(Tc = O.lSkms^^. A noise correction for CTc (Miesch & 



Bally|[l994| has negligible effect. The dispersion m cen 
troid values at a given lag is used as the error estimate 
in S* at that lag. 
Figure 20 shows S*{t) computed for C^^O over the 



high resolution region of Figure 10 There is a shallow 
increasing slope for t < 0.2 pc. A power-law fit to S*{t) 
in this region yields a slope of 0.24 ± 0.01. The function 
S'v(t), computed according to Equations [o] and 15 has 
amplitudes much greater than S'p(t) everywhere m our 
map implying that density fluctuations do not dominate 
on any scale. 

7.1.4. Discussion 

The values of a derived from the VGA analysis, the 
line width-size relationship, and the second-order veloc- 
ity structure function are, 0.72 ± 0.1, 0.14 ± 0.02, and 
0.24 ±0.01, respectively. These values differ signiflcantly 
according to the stated errors, but systematic errors may 
be important. 

The VGA analysis is sensitive to gas temperature 
(Lazarian, A. 2006, private communication), and there- 
fore varying kinetic temperatures throughout the cloud 
could confuse the VGA results. Gas temperatures are 
observed to vary between 9 and 25 K. Additionally, the 
slope of the two-dimensional power spectrum for a thin 
slice derived from Figure 18 is only a lower limit, since 



the trend toward larger slope indices may continue to 
smaller velocity widths than we have sampled with our 
data. However, a larger power spectrum index would 
only make the derived a values larger, and hence more 
discrepant. 

The distribution of S* (r) values at each r is skewed due 
to the full spread in values being of order the mean at 
most lags. Using the median instead of the mean steep- 
ens the slope considerably to a = 0.47 ± 0.02. However, 
this still does not reconcile the disparate derived values 
for a. 

Lastly, the consistency of these analyses depends upon 
the turbulence being homogeneous and isotropic. The 
structure seen in the ^^GO line wings in §[6]suggests that 
the energy input into the L1551 cloud is neither homoge- 
neous nor isotropic, and it is not clear if this is a sound 
assumption for the turbulence in the cloud. 

7.2. The Dissipation and Replenishment of Turbulent 

Energy 

The total turbulent energy in L1551 is a substantial 
fraction of the gravitational energy. However, turbulent 



energy dissipates on short timescales ( [Stone et al.||1998^ 
Mac Low fc Kless"eri||2004 1 suggesting that the observec 
turbulence in L1551 is not primordial. 

Galactic shear is negligible at the galactocentric dis- 
tance of L1551 given its size, and no evidence of large 
scale rotation is seen in the first moment map of ^"^GO 
emission. There are no signs of recent supernova explo- 
sions in the proximity of L1551, nor are there Hii re- 
gions nearby. The ^^GO line wings show coherent, thin 
structures that are expected a t the early time evo lution 

pOTter et al.|[l992 l. 



This 



of supersonic turbulence {e.g. 

structure is not random, but is clearly associated with 
the young stars in the region. 
The total energy in the "'^^GO line wings estimated in 



6.4 



is 4- 
and 



-10 X 
and 



1044, 

the 



ergs. Both the dyna mical t imescale 
diffusion timescale ( 
line wing features are ~ 10^ years, and 
L1551 IRS5 outflow has been estimated to be the same 



5.3.21 
tTle 



for the 
age of the 



order of magnitud e (Moriarty-Schieven & Snell 1988 



Richer et al. 2000 1 . i'herefore the observed energy m 



the cloud has been injected over the past ^ 10^ years, 
and the observed rate of energy input into LI 551 is 
i^input « 2 X 1032ergss-\ or 0.05 2.©. 

The energy from stellar feedback is expected to be 
highly turbulent and therefore the turbulent dissipation 
timescale, tdiss = £'turb/£'diss, sets the rate at which this 
energy is lost from the cloud. The dissipation of turbu- 
lent energ y depends on the scale at which that turbulence 



Ostrikerl 



is driven (Mac Low 1999 
|1996p. 



integrated 



igure 



5T 



1998 Gammie & 



Stone et al. 

shows the power spectra of the 
O intensity in the red and blue line wings. 
The overall best flt for a single power law is shown as the 
thin lines. The slopes of these spectra are more shallow 
than the slopes in the line core, and there are also clear 
deviations from the linear trend at scales ^ 1'. Since 
the structures that are clearly associated with embed- 
ded sources dominate the emission in the line wings, the 
broad excess peaks in the line wing power spectra signify 
a preferential scale at which energy is being injected into 



(arcmin) 




Fig. 21. — Intensity power spectra for the integrated intensity in 
the red (upper) and blue (lower) line wings of the ^''CO emission. 
The power spectrum of the blue line wing is divided by 100 for 
viewing purposes. The velocity intervals for the integrated inten- 
sity are 7.5-9.0 kms~^ and 3.4-5.7kms~^ for the red and blue line 
wing data, respectively. Clear excursions from the overall power- 
law form (lines) are seen near 1' scales. 
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the cloud. 

Using the dispersion of the composite C^^O spectrum, 
0.31 kms~^, we estimate the flow crossing time at an en- 
ergy injection scale of 0.05 pc, tf(Apoak) ~ 0.1 Myr. This 
timescale sets the overall sc aling for idiss- In the nu- 
merical studies by Stone et al., (1998j ) , tdiss was found to 
vary between ~ 0.5-0.8 1{ for a range in magnetic field 
strength between and ^ 50 /iG using a kinetic tem- 
perature of 15 K and a particle density equal to (n) w 
lO^cm-^. For Uiss = 0-7 U and E'turb = 4.3 x lO'^'ergs, 
the turbulent energy dissipation rate in L1551 is i?diss ~ 
2 X 10'^^ ergs s~^, approximately equal to the estimated 
energy input rat e from outflow. Using Equation 7 from 
Mac Low ( 1999 ) with the appropriate input values gives 



a consistent number for i?diss- 

The uncertainty in the stellar feedback energy is 
roughly a factor of 2 and the timescale over which the 
observed energy in outflow has been injected into the 
cloud is accurate to roughly the same level. The turbu- 
lent energy figure is relatively accurate, but our derived 
dissipation timescale may be off by a factor of 2 or more 
based on the broadness of the peaks seen in Figure [21] 
and the uncertain magnetic field strength in L1551 (see 

§|8]). Therefore we can only conclude that -Ediss ~ £'input 
to order of magnitude. However, this rough agreement 
between energy rates along with the lack of evidence for 
any other form of energy input, the age of LI 551 be- 
ing several dynamical times, and the fact that turbulent 
energy only survives for roughly a flow crossing time, 
strongly suggest that the observed turbulent energy has 
been supplied by outflows from forming stars. 

8. DISCUSSION 

The evolution of the L1551 cloud is influenced signif- 
icantly by the stars forming within it. Ample evidence 
for energy injection and cloud destruction is apparent in 
our high-resolution maps. The observed turbulent mo- 
tions in the cloud are short lived, ^diss ^ 0.1 Myr, and 
exist primarily on small scales where they may act as a 
turbulent pressure to s upport the cloud against gravity 
(Bonazzola et al. 19921. Additionally, the energy input 
rate is estimated to be of the order of the dissipation rate. 
These calculations all support the idea that the progres- 
sion of star formation is signiflcantl y influenced by the 
forming stars {e.g. Norman & Silk 1980 Li & Nakamura 
2006] ). 

The observed star formation efficiency, SFE, is de- 
fined as the present day stellar mass divided by the total 
present day mass in stars and gas. In L1551, SFE « 12%. 
The total star formation efficiency, SFEtot, is the fraction 
of total initial cloud mass that will ever become stars and 
may be as low as 9% in L1551 if the initial mass of the 
cloud were 1.6 times the mass of the present day cloud. 
It appears that L1551 is not capable of producing signif- 
icantly more stars in the future limiting SFEtot ^ 15%. 
However, further study is needed in the high column re- 
gions to the northwest. The theoretical expectation of 
SFEtot from o utflow-regulated star form ation is in the 



20001 



range 25-75% ( [Matzner fc McKee 

The age spread of PMS stars, fractional numbers of 
protostars, wTTSs, and cTTSs, and the spatial and 
kinematic distributions of the stellar component suggest 
LI 551 has been forming stars for much longer than its 



present day dynamical time. The near equality of stellar 
feedback and turbulent dissipation provides a means for 
the cloud to stay close to vir ial equilibrium over its stel- 



lar pr oduction history {e.g., Tan et al. 2006 Krumholz 



20061. However, the rate of star formation has not re- 



mamed constant in L1551. Roughly 9 protostars have 
formed in L1551 in the past ^ 0.1 Myr. If this rate of 
star formation were projected backward in time even a 
few million years, the number of stars produced in L1551 
would be in the hundreds. Therefore, we conclude that 
the increase in the rate of star formation over the life- 
time of L1551 is robust. Evidence for accelera ting star 
formation ha s been seen in Taurus as a whole (Palla & 
Stahler] 2002[) and in other w ell-studied star formmg re- 
gions (P& M fc Stahler|2000| [Huff fc Stah"Ier||2006[ ) and is 
expected it the star formati on rate follows the dyn amical 
rate of a contacting cloud ( Palla fc Stahler|[l999 l. 

Accelerating star formation fits loosely wTthin the 
frame work of dynamical star formation {e.g., Elmegreen 
20001. However, dynamical models typically neglect or 
discount the effects of stellar feedback. This is untenable 
in the case of L1551. The timescales of purely dynamical 
star formation, typically l-2tdyn, also do not match our 
observations. Even in the likely case that the dynamical 
timescale of the initial cloud was longer than the present 
day estimation, clear evidence for the formation of an- 
other star or stellar system that has yet to form pushes 
the productivity lifetime of L1551 well beyond the limits 
of purely dynamical formation. 

Missing from this picture of star formation is infor- 
mation regarding the magnetic field. Polarization mea- 
surements of background starlight in the periphery of 
L1551 reveal a magnetic field direction somewhat sug- 
gestively aligned with the L1551 IRS5 outflow direction. 
However, a sub- millimeter polarization measurement of 
the dusty clump around L1551 IRS5 suggests a fi eld di- 
rection nea rly perpendicular to the outflow axis (Vallee 
fc Bastien 2000). And while a firm lower limit of 2 //G 
has been set for the magnetic field strength in LI 551 
via rota tion measures of bac kground extragalactic radio 
sources (Simonetti fc CordeSjl986), useful constraints on 



the field strength are precluded without detailed infor- 
mation about the ionization state of the gas. 

The relative importance of the magnetic field in deter- 
mining the structure and evolution of L1551 is presently 
not clear. More detailed observations regarding the mag- 
netic field strength and orientation in L1551 are neces- 
sary to compliment this and other work {e.g., Moriarty- 



Schieven et al. 2006 Stojimirovic et al. 2006p and to 



stars {e.g., Shu et al. 


19871 IMcKee 


Mouschovias]|l993 


Shu et al. |2004|). 



1989 Ciolek fc 



9. CONCLUSION 

This article outlines an extensive observational case 
study of an active star-forming molecular cloud and 
its pre-main sequence association of stars. New, key 
observational additions to the study of this region in- 
clude wide-field, high-resolution maps of ^'^00(1 — 0) and 
C^^0(1 — 0) created through the combination of single- 
dish and interferometric data; deep, wide-field, near- 
infrared imaging; Spitzer IRAC data; and the culling of 
2MASS and literature data regarding the young stellar 
population in the region. These data shed light on the 
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relationship between the diverse phenomena observed in 
L1551 and the progression of star formation in this cloud. 

Within a 1?6 radius of the C^^O intensity weighted 
mean position of L1551, 4'^31'^24% +18°10'00" (J2000), 
35 PMS stars are selected to be part of the L1551 as- 
sociation. The distribution of spectral class and derived 
bolometric luminosities suggest that star formation be- 
gan > 6 Myr ago. The spatial and kinematic distribution 
of PMS stars are also consistent with a long star forma- 
tion history. The relatively high current rate of star for- 
mation is clear evidence that the star formation rate has 
increased over the productivity lifetime of L1551, consis- 
tent with accelerated star formation scenarios. The star 
formation appears to be progressing east to west with in- 
dividual sites of star formation likely occurring in regions 
of quiescent gas. 

The L1551 cloud contains 160 Af0 of total mass as 
traced by dust extinction, factors of 1.4 to 2 higher than 
previous estimates from molecular emission. The dis- 
crepancies are likely due to the saturation of CO emis- 
sion in high column regions and the destruction of CO 
in areas of low extinction. The dynamical timescale for 
the cloud, given the radial distribution of mass within 
its 0.9 ulius, is tdyn = 1.1 Myr. This is significantly 
shorter than the inferred age of L1551 based on its stellar 
component. 

The composite C-'^^O line width infers a total kinetic 
energy in the cloud in accord with virial equilibrium as- 
sumptions. The turbulent velocity field in L1551 pro- 
duces a density power spectrum Pp{k) oc k~^'^. We de- 
rive an energy spectrum index /3 « 1-2, but the observed 
turbulent motions exist primarily on small scales as in- 
ferred directly from the line width-size relationship in the 
cloud. The decay timescale for this turbulent energy is 
short, tdiss ^ 0.1 Myr and is not primordial in origin. 

The high-resolution ^^CO emission traces well the stel- 
lar feedback in the cloud. Thin, filamentary features in 



the line wings of the ^'^CO emission signify the interface 
between outflow and quiescent gas and are used to es- 
timate the stellar feedback, i^input ~ 0.05 L©. This is 
not the full mechanical luminosity of outflow in L1551, 
rather a signiflcant fraction of outflow energy contributes 
to the destruction of the cloud or is otherwise lost to the 
greater interstellar medium. It is estimated that 4-5 Mq 
of gas has been excavated from the cloud by the current 
generation of outflows. The preferred scale at which en- 
ergy is injected into the cloud is Apeak ~ 0.05 pc, about 
a thirtieth of the diameter of the cloud as seen in dust 
extinction. 

The dissipation rate of turbulence and the rate of en- 
ergy injection from outflow balance according to our cal- 
culations. There are no other significant sources of en- 
ergy seen in L1551, and therefore outflows provide a sig- 
nificant supply of kinetic energy and play an important 
role in the evolution of the cloud. 
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Properties of Selected Pre-main Sequence Stars 
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